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^ 1 Introduction 

(D , 
1-^ . Two of the most challenging tasks in molecular simulation consist in capturing the 

properties of systems with long-range interactions (e.g. electrolyte solutions), and of 
systems containing large molecules such as hydrogels. These tasks become particu- 
larly demanding when explicit solvent models are used. Therefore, massively parallel 

^.,^. supercomputers are needed for both tasks. 

(~| ' For the development and optimization of molecular force fields and models, a large 

Q^, number of simulation runs have to be evaluated to obtain the sensitivity of thermody- 

namic properties with respect to the model parameters. This requires both an efficietn 
work flow and, obviously, even more computational resources. The present work dis- 

^ ' cusses the force field development for electrolytes regarding thermodynamic prop- 

^^O , erties of their solutions. Furthermore, simulation results for the volume transition of 

hydrogels in solution containing electrolytes are presented. Both applications are of 
interest for engineering. It is shown that the properties of these complex systems can 
be reasonably predicted by molecular simulation. 



2 Development of force fields for alkali and halogen ions in 
aqueous solution 

2.1 Outline 

The simulation of electrolyte systems is computationally very expensive due to the 
long-range interactions, which have to be taken into account by suitable algorithms. 
Examples of such algorithms are the classical Ewald summation [1] and its deriva- 
tives like particle mesh Ewald summation [2] or the particle particle/particle mesh 
method [3]. 
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Early efforts in this research area were mainly directed to the development of ion 
force fields, which were capable of reproducing static structural properties of solu- 
tions [4, 5, 6]. Recently, however, these models have been proven to be too inaccurate 
for the prediction of basic thermodynamic properties [7, 8]. Since the models were 
parameterized solely at short distances, long distance effects are underestimated. 
These effects can be seen in particular in the density and the activity of electrolyte 
solutions. 

In the present work new atomistic models for alkali and halogen ions are presented, 
which accurately describe not only structural properties of aqueous electrolyte solu- 
tions but also basic thermodynamic properties like their density. The main focus of 
the model development is placed on transferability, i.e. the ion models are intended 
to describe solution properties independently of the cation/anion combination. 
The ions were modeled as Lennard-Jones (LJ) spheres with superimposed charges 
of +1 in units of elementary charges, located in the center of mass. For the solvent 
water, the SPC/E model [9] was used, which consists of one LJ site and three point 
charges. 



2.2 Simulation details 

All simulations in Section 2 were performed using the Monte-Carlo technique in 
the isobaric-isothermal (NpT) ensemble at 20 °C and 0.1 MPa. The simulation vol- 
ume contained a total of 1000 molecules and ions, respectively. The electrostatic 
long-range interactions were calculated using the Ewald summation with an Ewald 
parameter K of 5.6/L, where L denotes the length of the cubic simulation volume. 
The dispersive and repulsive long-range contributions were approximated using the 
assumption of an homogeneous fluid beyond the cut-off distance of 1 1.9 A. The sim- 
ulation program employed was an extended version of ms2 [10]. 



2.3 Solvent model and reduced properties 

For the calculation of aqueous electrolyte solutions, the model for water is crucial, 
since it represents the largest fraction in the mixture. The SPC/E [9] water model 
was employed, which reproduces the density and other thermodynamic properties at 
chosen conditions in good agreement with experimental data (Pspc/e{T — 2Q°C,p — 
O.lMPa) = 999.5 g/1). To minimize the influence of errors in the solvent model on 
the parameterization of the ion force fields, the reduced density p of the aqueous 
solution was chosen as objective function, which is defined as the fraction of the 
density of the electrolyte solution Pes and the density of the pure solvent ps at the 
same temperature and pressure 

P = ^. (1) 

Ps 

Figure 1 shows a typical plot of the reduced density of an aqueous electrolyte so- 
lution as a function of salt mass fraction using sodium chloride as an example. The 
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dependence of p with x™ is almost linear. It can also be seen from Figure 1 that elec- 
trolyte force fields from the literature predict this linear coiTelation, but with a wrong 
slope. Note also that the solvent activity (or activity coefficient) can be regarded as 
a normalized property. However, a similar approach fails for properties like the ion 
self-diffusion coefficient. 



"^ 




0.20 



Fig. 1. Reduced density p as a function of sodium chloride mass fraction. The line indicates 
the experimental values [11], while the symbols represent simulation results using varying ion 
force fields [12, 13, 14, 15]. 



2.4 Ion force field 



A global parameterization strategy was used in the present study. Due to the simple 
model for the ions, the overall parameter space is small. It contains the two LJ param- 
eters a and e for each ion. The position and magnitude of the charges are constant. 
In a preliminary analysis, the sensitivity of the reduced density on the Lennard- Jones 
parameters was evaluated. It shows that the influence of the LJ energy parameter e 
on p is very small. Therefore, the a parameters were determined by a fit to the re- 
duced density, while the £ parameters were set here to a constant value, i.e. 200 K 
for sodium and chloride, based on results of a study on the water activity in the elec- 
trolyte solution, which is not described in detail here. 

The objective function of the optimization used to determine the LJ cr parameters is 
the slope of the reduced density p with increasing salt mass fraction at x'™' -^ 



dp 

djcf™: 



a:("i)^0 



( dp ^ 



(2) 



4 J. Walter, S. Deublein, S. Reiser, M. Horsch, J. Vrabec and H. Hasse 

In molecular simulations, these slopes were systematically derived for varying values 
of aanion ^nd Ocation- The data were smoothed by a two dimensional polynomial fit. 
The ion parameters for all alkali and halogen ions were chosen such that the sum 
of the squared deviations between the functional fit and the experimental data is 
minimized. 

2.5 Results 

New atomistic models for alkali and halogen ions were determined that describe 
the reduced density of the aqueous electrolyte solution in good agreement with ex- 
periments. As an example the force field for sodium chloride is shown in Table 1. 
Figure 1 shows the excellent agreement between the simulation using the new model 
and the experimental results. 



Table 1. Force field of sodium chloride 



Ion CT / A £ / K 

Na+ 1.88 200 

C\- 4.41 200 



The new atomistic models for alkali and halogen ions were investigated regarding 
the representation of structural properties of aqueous solution. Radial distribution 
functions gij(r) were used for the characterization which are defined by 

gijin ^ - — ■ (3) 

Pj.bulk 

Here pj (r) is the density of component j as a function of the distance r between two 
ions of component / and j, respectively, and p^.buik is the number density of compo- 
nent j in the bulk phase. For the characterization of aqueous electrolyte solutions, 
the radial distribution function of water around the ions is of particular interest. In 
this case, water is represented by the position of the oxygen atom. For a solution of 
sodium chloride in water, ^,.h2o('') is shown in Figure 2. 

The locations of the first maximum and minimum, respectively, for both ions are 
in good agreement with experimental measurements [16], as shown in Table 2. The 
hydration shell of the cation is more attached to the ion than of the anion, which is 
represented by the height of the first peak in the distribution function. 
The attraction of the solvent to the cation is also visible in the hydration number 
around the ions, which describes the number of solvent molecules in the closest 
vicinity to the solute / and is defined by 

Wi.HzO ^ 4;rpH20 / '" ^i,H2o('-)dr . (4) 
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Fig. 2. Radial distribution function gi,H2o('') of ths solvent water around the sodium cation 
and chloride anion, respectively, at T =20 °C and p =0.1 MPa. 



Table 2. Location of the first maximum rmax and minimum r^^^^ of the radial distribution func- 
tion g,;H,o('') for an aqueous sodium chloride solution. The simulative results are compared 
to experimental data from the literature [16]. 



Ion / 



Exp , 
' max ' J 



„Sim / 



r^l' 



Na^ 



2.2 
3.4 



2.3 
3.3 



3.0 
3.9 



3.0 
4.0 



Here, p, defines the number density of component / and rmin the distance up to which 
the hydration number is calculated. For the first shell, the value of rmin is chosen to 
be the distance of the first minimum in the radial distribution function. For sodium 
chloride for example, the calculation of hydration numbers reproduces experimental 
results in good agreement, cf. Table 3. 



2.6 Computational demands 

Typical simulations to generate data points for the studies discussed in Section 2 were 
carried out on 16 CPUs running for 72 hours. For the prediction of other thermody- 
namic data like activity coefficients, up to 32 CPUs running for 72 hours depending 
on the system are required. For these simulations a virtual ram of 216 MB was re- 
quired. 
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Table 3. Hydration number n of Na+ and Cl^ in aqueous solution. The simulative results are 
compared to experimental data from the literature [16]. 



Ion 



Na+ 5.6 5-6 

C\- 7.5 7-8 



3 Self-diffusion coefficients of solutes in electrolyte systems 

3.1 Outline 

The self-diffusion coefficient is, in comparison to the density, an individual property 
of the different solute species and the solvent in electrolyte systems. Self-diffusion 
coefficients of anions and cations in aqueous solution are experimentally accessible 
and numerous experimental data are available in the literature [17, 18]. Therefore, it 
is worthwhile trying to fit model parameters of ion force fields to self-diffusion coef- 
ficient data. As a suitable strategy for determining the size parameters a is available 
(cf. Section 2), the question may be raised whether self-diffusion coefficient data are 
useful for determining the energy parameters e. 

In molecular simulations, the self-diffusion coefficient is usually determined by time 
and memory consuming methods, which require simulations of large systems. Ex- 
amples for such methods in equilibrium molecular dynamics are the mean square 
displacement [19] and the Green-Kubo formalism [20, 21]. In this formalism, the 
self-diffusion coefficient is related to the time integral of the velocity auto-correlation 
function. The calculation of self-diffusion coefficients of solutes in electrolyte sys- 
tems are computationally expensive due to additional time consuming algorithms 
(e.g. Ewald summation [1]) that allow for a truncation of ionic interactions in molec- 
ular simulations. 

In aqueous solution, the cations and anions are surrounded by a shell of strongly 
bonded water molecules (hydration shell). These hydrated ions diffuse within a bulk 
fluid, which is itself also highly structured. Therefore, the mobility and accordingly 
the self-diffusion coefficient of ions is strongly related to the structure of the water 
molecules around the ions [22]. 

3.2 Methods 

The investigated solution consisted of sodium and chloride ions as solutes and ex- 
plicit water as solvent. The ions were modeled as Lennard-Jones spheres with a cen- 
tral charge. Size parameters a for different ion force fields were determined as dis- 
cussed in Section 2. The energy parameter e was modified in a range of 50 to 250 K. 
Water models were taken from the literature. 
First, the density of the solution was determined in an isobaric-isothermal (NpT) 
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molecular dynamics (MD) simulation at a desired temperature and pressure. Then, 
the velocity auto-correlation function and, according to the Green-Kubo formalism 
[20, 21], the self-diffusion coefficients were determined in an isochoric-isothermal 
(NVT) MD simulation at this temperature and density. The sampling length of the ve- 
locity auto-correlation function was set to 1 1 ps and the time span between the origin 
of two auto-correlation functions was 0.06 ps. The separation between the time ori- 
gins was chosen such that all autocorrelation functions have decayed at least to 1/e 
of their normalized value. For both NpT and NVT simulations, the molecular dy- 
namics unit cell with periodic boundary condition included 4420 water molecules, 
40 sodium and 40 chloride ions. The long-range charge-charge interactions were 
calculated using Ewald summation [1]. The simulation program employed was an 
extended version of ms2 [10], which is developed by our group. 

3.3 Results 

The self-diffusion coefficients of anions and cations determined in molecular simu- 
lation largely depend on the used molecular model of water, as the mobility of the 
ions is influenced by the hydration shell and the structure of the surrounding water. 

Self-diffusion coefficient of water models 

The accuracy of the estimated self-diffusion coefficients of pure water for different 
water models is verified with respect to experimental data [23]. The self-diffusion 
coefficients were determined with the Green-Kubo formalism [20, 21] in molecular 
dynamics simulations at a temperature of 25 °C and a pressure of 0.1 MPa. For this 
study, three commonly used rigid nonpolarizable molecular water models of united- 
atom type, namely SPC/E [9], TIP4P [24] and TIP4P/2005 [25], were chosen. 



Table 4. Self-diffusion coefficients of pure water at 25 °C and 0. 1 MPa of different molecular 
water models. The number in parenthesis indicates the uncertanty of the last digit. 

Model SPC/E TIP4P TIP4P/2005 Experiment 

L»w [m^s-i] 26.2(1) 36.7(2) 21.9(1) 22.3 (-) 



The determined self-diffusion coefficients of the different water models vary over a 
wide range, cf. Table 4. For TIP4P/2005, the self-diffusion coefficient of pure wa- 
ter is in good agreement with the experimental value. In contrast, both the SPC/E 
and the TIP4P model overestimate the mobility of water molecules in pure water. 
The obtained values for the self-diffusion coefficient are in good agreement with the 
results published by other authors [26] . 
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Model development 

For fitting the energy parameter e of the ion force fields to experimental self-diffusion 
coefficients of anions and cations in aqueous solution, the influence of e on the deter- 
mined self-diffusion coefficients in molecular simulations was investigated. In this 
study, the above mentioned three water models were used. 

However, it turned out that the number for the energetic parameter £ has no signif- 
icant influence on the self-diffusion coefficient. This study is not discussed here in 
detail. In the following, e = 200 K is used. The results for both sodium and chloride 
are shown in Figure 3. 
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Fig. 3. Self-diffusion coefficient of sodium (Z^Na) ^^'^ cliloride (Dq^) in aqueous solution at 
25 °C and 0.1 MPa for different water models. The energy parameter e for both sodium and 
chloride force fields was set to 200 K. Experimental data for sodium [17] and chloride [18] 
were taken from the literature. 



As can be seen in Table 4 and Figure 3, there is no correlation between the accuracy 
of the determined self-diffusion coefficient of pure water and the accuracy of the esti- 
mation of ion mobility for the different water models. For example, the TIP4P model 
significantly overestimates the self-diffusion coefficient of pure water, whereas the 
simulation results for the self-diffusion coefficients of the sodium and the chloride 
ion are in fair agreement with experimental data (deviation for sodium of 4%; devi- 
ation for chloride of 18%). 

All water models overestimate the interaction between the water molecules in the 
hydration shell, the sodium ion and the bulk fluid. Hence, the cation mobility is too 
low. The same is true for the chloride ion only for SPC/E and TIP4P/2005 water 
models, whereas the TIP4P underestimates the interaction in the hydration shell. 
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3.4 Computational demands 

All molecular simulations in Section 3 were carried out with the MPI based molec- 
ular simulation program ms2, which is developed in our group. The total computing 
time for determining self-diffusion coefficients of ions in electrolyte systems was 
216 hours on 36 CPUs (72 hours for the NpT run and 144 hours for the NVT run). 
These simulations require large systems as the accuracy of the Green-Kubo formal- 
ism for ions increases with increasing number of solutes and, at the same time, an 
infinite dilution is aspired. For these simulation a maximum virtual memory of 1 .76 
GB was used. 



4 Hydrogels in electrolyte solutions 
4.1 Outline 

Hydrogels are three-dimensional hydrophilic polymer networks. Their most char- 
acteristic property is their swelling in aqueous solutions by absorbing the solvent, 
which is influenced by various factors. Hydrogels can be used in many applications 
like e.g. superabsorbers such as in diapers [27] and contact lenses [28]. To fully 
exploit the potential of hydrogels in all these applications, it is important to under- 
stand, describe and predict their swelling behavior The hydrogel which is studied 
in the present work is built up of poly(N-isopropylacrylamide) (PNIPAAm) cross- 
linked with N,N'-methylenebisacrylamide (MBA). PNIPAAm is one of the most 
extensively studied hydrogels in the scientific literature and is mainly used in bio- 
engineering applications [29] . The degree of swelling in equilibrium of PNIPAAm 
is significantly influenced by many factors [30, 31, 32, 33]. On the one hand, the 
swelling depends on the structure of the hydrogel itself, like the type of the monomer, 
but also the amount and type of cross-linker and of co-monomers. On the other hand, 
the environment conditions like temperature, type of solvent, solvent composition, 
electrolyte concentration or pH-value of the solvent influence the swelling. Varying 
these factors, the hydrogel typically shows a region where it is swollen and a region 
where it is collapsed. In between those two regions lies the region of volume tran- 
sition. The solvent composition which is characteristic for that transition is called 
©-solvent here. The ©-solvent mainly depends on the environmental factors and the 
nature of the polymer chain but not on the amount of cross-linker. 
For the quantitative description of the swelling of hydrogels, various types of models 
are used [34]. It is normally not possible to quantitatively predict the swelling of hy- 
drogels or its dependence on factors the models were not adjusted to. With molecular 
simulation it is possible to predict the swelling of different hydrogels upon varying 
any environmental factor, as was shown in a previous study [35]. 
In the present work, the swelling of PNIPAAm hydrogel is studied with atomistic 
molecular dynamics simulation. The results for the ©-solvent are compared to ex- 
perimental data of PNIPAAm hydrogel as a function of the temperature in water 



10 J. Walter, S. Deublein, S. Reiser, M. Horsch, J. Vrabec and H. Hasse 

[36]. The experimental data shows that the ©-solvent of PNIPAAm in different elec- 
trolyte solutions follows the Hofmeister series. For this study the electrolyte sodium 
chloride (NaCl) was considered. 

4.2 Models 

For the molecular dynamics simulations of PNIPAAm in aqueous solutions, the 
OPLS-AA force field [37, 38] was employed to describe PNIPAAm. It was com- 
bined with the SPC/E water model [9]. In previous studies, it was shown that this 
combination allows predicting the volume transition of PNIPAAm in water as func- 
tion of the temperature [35]. Different NaCl models from the literature were used 
and compared to the NaCl model developed in this work. The used electrolyte mod- 
els from the literature are GROMOS-96 53 A6 (G53A6) [39] and KBFF [40]. 



4.3 Simulation details 

Molecular simulations of PNIPAAm single chains were carried out with version 4.0.5 
of the GROMACS simulation package [41, 42]. Simulations with PNIPAAm in aque- 
ous NaCl solutions at 25 °C were performed in order to find the best NaCl model for 
the simulation of the volume transitions of PNIPAAm in NaCl solutions. 
In a previous study, it was shown that the amount of cross-linker of the hydrogel 
has no significant influence on the ©-solvent of the hydrogel and that this value is 
the same as for the PNIPAAm polymer [35]. Therefore, the simulations can be per- 
formed with a single PNIPAAm chain. For equilibration, single PNIPAAm chains 
in water were simulated in the isobaric-isothermal ensemble (NpT) over 1 to 5 • 10^ 
timesteps. The pressure was 0.1 MPa and was controlled by the Berendsen barostat 
[43], the temperature was controlled by the velocity rescaling thermostat [44] and the 
timestep was 1 fs for all simulations. Newton's equations of motion were numerically 
solved with the leap frog integrator [45]. For the long-range electrostatic interactions, 
particle mesh Ewald [46] with a grid spacing of 1.2 A and an interpolation order of 
four was used. A cutoff radius of re = 15 A was assumed for all interactions. After 
equilibration, 2 to 6 • 10^ production time steps were carried out with constant simu- 
lation parameters. Note that the production steps include the conformation transition 
as well as the simulation of the equilibrium. 
In order to analyze the results, the radius of gyration Rg was calculated 

V ^if^i ) 

which characterizes the degree of stretching of the single chain, where m, is the mass 
of site / and ||r,|| is the norm of the vector from site / to the center of mass of the 
single chain. The radius of gyration in equilibrium was calculated as the arithmetic 
mean over the last 5-10^ time steps of the simulation together with its standard 
deviation. 
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4.4 Results and discussion 

For the temperature of 25 °C, the hydrogel is swollen in pure water and collapsed 
above the NaCl concentration x{^q of about 0.03 g-g^' [36]. Therefore, the sin- 
gle chain should be stretched in pure water and low electrolyte concentrations and 
collapsed at electrolyte concentrations above 0.03 g-g^^. Figure 4 shows the simula- 
tion results as the radius of gyration over the NaCl concentration. The NaCl model 
from this work and G53A6 are able to predict the ©-solvent of PNIPAAm at a NaCl 
concentration of 0.03 gg^'. The NaCl model KBFF does not yield the fully col- 
lapsed conformation. The models, that show the ©-solvent at 0.03 gg^' also show 
a more or less stretched single chain at electrolyte concentrations above this point. 
The model with the best results here is the one from this work, for it only yields a 
slightly stretched conformation of the chain at electrolyte concentrations above the 
©-solvent. 
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Fig. 4. Radius of gyration Rg of a PNIPAAm chain of 30 monomers in NaCl solutions of 

about 14,000 solvent molecules in equilibrium as a function of the NaCl concentration x^L 
for the different NaCl models at a temperature of 25 °C. The error bars indicate the standard 
deviation. 



It is therefore clearly possible to obtain qualitative predictions for the swelUng of 
PNIPAAm hydrogels in electrolyte solutions of NaCl by molecular simulation of a 
single chain. For the ©-solvent in NaCl solutions, it was even possible to quantita- 
tively reproduce experimental data. First studies on the volume transition of PNI- 
PAAm hydrogels in electrolyte solutions of sodium sulfate Na2S04 were currently 
carried out (results not shown here). By comparing the results of the two electrolyte 
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solutions NaCl and Na2S04 in water, it is also possible to determine the effect of the 
Hofmeister series on the solubility of PNIPAAm in aqueous electrolyte solutions. 
The Hofmeister series leads to a ©-solvent of PNIPAAm in Na2S04 solutions at a 
lower concentration than in NaCl solutions [36]. This correlation could be repro- 
duced by the molecular simulations. In summary, this is an unexpectedly favorable 
agreement between predictions by molecular simulation and experimental data, es- 
pecially when considering that the force fields were not adjusted to any such data. 



4.5 Computational demands 

All simulations presented in Section 4 were carried out with the MPI based molec- 
ular simulation program GROMACS. The parallelization of the molecular dynamics 
part of GROMACS is based on the eighth shell domain decomposition method [42]. 
With GROMACS, typical simulation runs to determine the radius of gyration in equi- 
librium employ 128 CPUs running for 24-72 hours. For these simulations very large 
systems must be considered comprising typically about 58 800 interaction sites. For 
these simulations a maximum memory of 284 MB and a maximum virtual memory 
of 739 MB was used. 



5 Conclusion 

This work covers the development of ion force fields for describing the thermody- 
namic properties of electrolyte solutions and the applications of these force fields 
for predicting the volume transition of hydrogels by atomistic molecular simulations 
with explicit solvent models 

The present work proves that alkali and halogen ions can be reliably modeled by a 
LJ sphere with a superimposed charge located at the center of mass. The developed 
force fields allow predicting structural properties like the radial distribution function 
and the hydration number as well as thermodynamic properties like the density. 
The self-diffusion coefficient of pure water was determined using various well known 
molecular models of water The agreement with experimental data is often poor Fur- 
thermore, ion self-diffusion coefficients were determined in simulations using these 
different water models. No correlation was observed between the accuracy of the 
self-diffusion coefficients of pure water and the accuracy of the self-diffusion coeffi- 
cients of the ions. Further research effort in developing a new water model is needed 
to gain accurate preditions of transport properties in aqueous electrolyte systems. In 
addition, the study shows that the self-diffusion coefficient of ions in aqueous solu- 
tions is almost independent of the LJ energy parameter e of the ion. 
With the developed electrolyte models, it was possible to predict the volume transi- 
tion of hydrogels in electrolyte solutions qualitatively and in some cases even quanti- 
tatively. The results also reproduce the effect of the Hofmeister series on the swelling 
of hydrogels. 
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